home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / KSMIRNOV.C < prev    next >
C/C++ Source or Header  |  1996-03-31  |  3KB  |  119 lines
  1. /* ============ */
  2. /* ksmirnov.c    */
  3. /* ============ */
  4. #include <math.h>
  5. #include <mconf.h>
  6. extern    double    MINLOG;
  7. /* ==================================================================== */
  8. /* KSmirnov - Calculates Kolmogorov-Smirnov Probability for Given n & e */
  9. /* ==================================================================== */
  10. double
  11. KSmirnov(int n, double e)
  12. {
  13.     int     v, nn;
  14.     double  evn, omevn, p, t, c, lgamnp1, RetVal;
  15.  
  16. /* ------------------------------------------------------------ *
  17.  * Exact Smirnov statistic, for one-sided test [see ref].    *
  18.  *                                *
  19.  *            [n(1-e)]                    *
  20.  *    +           -        v-1         n-v    *
  21.  *  Pr{D <= e} = 1 - e >  C   (e + v/n)    (1 - e - v/n)    *
  22.  *    n           - n v                    *
  23.  *              v=0                    *
  24.  *                                *
  25.  *  This is the probability that the statistic Dn <= e.     *
  26.  *                                *
  27.  *  [n(1-e)] is the largest integer not exceeding n(1-e).    *
  28.  *                                *
  29.  *  nCv is the number of combinations of n things taken v    *
  30.  *    at a time (binomial coefficient).            *
  31.  *                                *
  32.  *  Reference:                            *
  33.  *  Z.W. Birnbaun and Fred H. Tingey, One-Sided Confidence    *
  34.  *  Countours for Probability Distribution Functions, Annals of *
  35.  *  Mathematical Statistics 22(1951), pp 592-596        *
  36.  *                                *
  37.  *  Coded by: Stephen L. Moshier, October 1995            *
  38.  *  Modified by: K. B. Williams, December 1995            *
  39.  * ------------------------------------------------------------ */
  40.  
  41.     if (n <= 0 || e <= 0.0 || e > 1.0)
  42.     {
  43.     if (e == 0)
  44.     {
  45.         RetVal = 0;
  46.     }
  47.     else
  48.     {
  49.         char    *ErrPrt;
  50.  
  51.         if (n <= 0)
  52.         {
  53.         ErrPrt = "KSmirnov (n <= 0)";
  54.         }
  55.         else if (e < 0)
  56.         {
  57.         ErrPrt = "KSmirnov (e < 0)";
  58.         }
  59.         else
  60.         {
  61.         ErrPrt = "KSmirnov (e > 1)";
  62.         }
  63.  
  64.         mtherr(ErrPrt, DOMAIN);
  65.         RetVal = -1;
  66.     }
  67.     }
  68.     else
  69.     {
  70.     nn = (int) floor((double) n * (1.0 - e));
  71.     p = 0.0;
  72.     if (n < 1013)
  73.     {
  74.         c = 1.0;
  75.         for (v = 0; v <= nn; v++)
  76.         {
  77.         evn = e + ((double) v) / n;
  78.  
  79.         if (evn < 1)
  80.         {
  81.             t = (double) (v-1) * log(evn) +
  82.             (double) (n-v) * log(1.0-evn);
  83.  
  84.             if (t > MINLOG)
  85.             {
  86.             p += c * exp(t);
  87.             }
  88.         }
  89.  
  90.         /* Next combinatorial term; worst case error = 4e-15.  */
  91.         c *= ((double) (n - v)) / (v + 1);
  92.         }
  93.     }
  94.     else
  95.     {
  96.         lgamnp1 = lgam((double) (n + 1));
  97.         for (v = 0; v <= nn; v++)
  98.         {
  99.         evn = e + ((double) v) / n;
  100.         omevn = 1.0 - evn;
  101.         if (fabs(omevn) > 0.0)
  102.         {
  103.             t = lgamnp1
  104.               - lgam((double) (v + 1))
  105.               - lgam((double) (n - v + 1))
  106.               + (v - 1) * log(evn)
  107.               + (n - v) * log(omevn);
  108.             if (t > MINLOG)
  109.             p += exp(t);
  110.         }
  111.         }
  112.     }
  113.  
  114.     RetVal = 1 - p * e;
  115.     }
  116.  
  117.     return (RetVal);
  118. }
  119.